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ABSTRACT 


A  numerical  scheme  which  has  been  successfully  employed  for 
two-dimensional  hydrodynamics  is  extended  to  the  axisymmetric 
visco-plastic  equations  governing  the  hypervelocity  impact  phenome¬ 
non,  Computing  procedures  are  outlined  which  are  complete  enough 
for  interior  regions  of  the  medium.  The  formulation  must  be  further 
developed,  however,  to  allow  for  such  contingencies  as  free  bound¬ 
aries,  interior  empty  cells  and  the  axis  of  symmetry  before  a  com¬ 
puting  code  can  be  written  for  the  impact  problem. 
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INTRODUCTION 


The  object  of  this  theoretical  investigation  of  the  hypervelocity 
impact  phenomenon  is  the  solution  of  the  governing  equations  for 
the  axisymmetric  case  depicted  in  Figure  1.  The  pertinent  visco¬ 
plastic  equations  were  presented  in  an  earlier  report  (Ref,  1). 
Numerical  schemes  have  been  developed  and  solutions  obtained  for 
th'e  one -dimensional,  time  dependent,  formulation  of  the  equations 
(Refs,  Z,  3);  the  results  have  been  amplified  in  other  discussions 
(Refs.  4,  5).  The  axisymmetric  problem,  however,  requires  the 
solution  of  the  two-dimensional,  time  dependent  equations,  A  large 
difference  exists  in  the  difficulty  between  the  1 -D  and  2-D  prob¬ 
lems. 


In  developing  a  finite  difference  formulation  of  the  axisymmet¬ 
ric  impact  problem,  extension  of  an  existing  scheme  devised  for 
two-dimensional  hydrodynamics  can  be  attempted  for  expediency. 
Several  methods  of  treatment  have  been  used  for  those  problems 
dependent  upon  two  or  more  space  coordinates.  These  variations 
usually  employ  (a)  I,agrangian  coordinates  in  which  the  mesh  of 
cells  is  imbedded  in  the  medium  and  moves  with  it  (Refs,  6,  7,  8), 
(b)  Euler ian  coordinates  which  are  not  fixed  in  the  medium  but 
are  usually  stationary  in  the  laboratory  frame  of  reference,  or  (c) 
a  mixed  Euler -Lagrange  system  which  attempts  to  take  advantage 
of  the  better  features  of  both  fixed  and  moveable  coordinates 
(Refs,  9,  10), 

The  difficulties  with  schemes  employing  Lagrangian  coordin¬ 
ates  is  the  large  distortions  which  are  involved  in  the  present 
problem.  The  Eulerian  systems  have  the  disadvantage  that  to 
account  for  the  projectile -target  interface  and  the  free  surfaces 
of  the  projectile  and  target  is  extremely  difficult.  These  are  the 
basic  reasons  for  the  decision  to  adopt  the  particle -in -cell  method 
(abbreviated  PIC)  which  has  been  developed  at  Los  Alamos  by 
Dr,  Francis  H,  Harlow  and  his  colleagues. 

The  essential  features  of  the  method  have  been  described  in 
two  Los  Alamos  reports  (Refs.  9,  10),  Briefly,  the  space  oc¬ 
cupied  by  the  medium  is  divided  into  a  suitably  chosen  mesh  of 
fixed  cells  through  which  the  medium  moves.  The  medium  within 
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each  cell  is  represented  by  a  set  of  mass  points  (particles).  At 
the  end  of  the  n  -th  time  cycle  the  mass  (equal  to  the  sum  of  the 
masses  of  the  particles  located  in  that  cell)  velocity,  pressure,  and 
specific  internal  energy  are  associated  with  each  cell.  To  obtain 
the  corresponding  data  at  the  end  of  the  (n  +  l)-st  time  cycle  one 
makes  a  three  phase  calculation.  In  Phase  I  the  field  functions  are 
changed  neglecting  the  motion  of  the  medium.  In  Phase  II  the  mass 
points  are  moved  and  the  field  functions  then  recalculated  to  account 
for  the  motion.  In  Phase  III  various  functionals  are  computed  which 
furnish  checks  on  the  accuracy  of  the  calculations, 

THE  PIC  FORMULATION 

The  projectile  is  considered  to  be  a  circular  cylinder  of  length 
L  and  diameter  H.  The  space  net  size  is  5r  by  Sz  as  depicted  in 
Figure  2,  For  example,  the  projectile  is  divided  into  £xd  cells, 

(1)  £  =  L/5z  d  =  H/2Sr 

each  of  which  is  actually  a  toroid  of  revolution.  The  cells  are 
labeled  with  index  ,  with  i  and  j  increasing  in  the  r  and  z 

directions,  respectively: 

(2)  /  j  X  /  (i  -  l)5r^  r<i  5r  ^ 
cell(^  i  /  ”  (  (j  —  1)  Sz  ^  z  <j  5z  / 

The  pressure  for  cell  is  pj  ;  the  average  pressure  along  the 

boundary  between  cells  and^.  |j^^is  denoted  by  P;  +  i/2’ 

the  average  pressure  between  cells  and  ^  is  p.  j 

Analogous  notations  are  used  for  the  other  cell-wise  quantities. 

Having  a  different  mass  for  different  mass  particles  is  con¬ 
venient  in  cylindrical  coordinates.  Each  particle  is  given  a  fixed 
mass  whose  value  is  proportional  to  the  radius  of  the  cell  within 
which  it  lies  originally  (t  =  0)  .  If  there  are  originally  k 
particles  in  each  cell,  then  the  mass  of  each  of  the  k  particles 
originally  in  cell  is 
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(3) 


2  ff  rj  St,Sz  2  fr  (i  -  1/2)  (5r)^  5z  p^  • 
k  k 

The  k  particles  are  assumed  to  be  originally  randomly  distributed 
in  the  cell.  The  location  of  each  particle  at  the  end  of  the  n  -  th 
time  cycle  is  stored  in  the  machine  memory  until  its  location  at 
the  end  of  the  (n  +  1)  -st  cycle  is  computed.  The  conservation  of 
mass  is  therefore  automatic  in  the  PIC  method, 

(a)  Phase  I  Calculations 


The  particles  are  not  moved;  thus  the  transport  terms  are 
dropped  from  the  momenta  and  energy  equations  governing  the 
visco-plastic  model  (see  Appendix  A),  The  tentative  new  cellwise 
velocity  components  and  specific  internal  energy  become: 
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where 
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(11) 


(r^)?  =  r  2 


2  u  r 
"o  o 


\/(D2)| 


(D2)i 


At  the  end  of  the  n  -  th  cycle  there  are  stored  in  the  memory 
the  mass,  specific  internal  energy,  velocity  components,  and  the 
pressure  for  every  cell.  Table  I  shows  these  together  with  the 
quantities  which  replace  them  during  the  sequence  of  computational 
steps  required  for  Phase  I  of  the  (n  +  1)  -  st  time  cycle  calculations; 


Step  1 ,  1  The  velocity  components  at  the  end  of  the  n  -  th  time 
cycle  are  evaluated  at  the  cell  boundaries: 


(12) 
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i 

i  -  1/2 


ul  .  +  uj 
1  —  1  1 


2 


j  -  1/2 


I  1 


2 


i  -  I  2  “  i 

(13)  2  2 

Step  1 , 2  The  divergence  of  the  velocity  at  the  cell  center  is 
computed  using  equation  (7), 

Step  1 . 3  The  value  at  the  cell  center  of  the  distortional 
strain-rate  invariant  is  computed  using  equation  (8), 

Step  1.4  The  value  at  the  cell  center  of  the  strain-rate  depen¬ 
dent  viscosity  coefficient  is  computed  using  equation  (9). 

Step  1 . 5  The  cell  center  values  of  the  quantities  P  and  Q  are 
computed  using  equations  (10), 
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Step  1,  6  Equation  (11)  is  used  to  compute  the  cell  center  value 
of  the  flow  statistic.  Further,  the  boundary  values  of  P,  Q  and  fi 
are  computed  according  to 

P|  .  +  pj  -  ^  +  p  j 

'  ~  ‘  -wo  '  i 


(14) 
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2 

Step  1,7  The  tentative  new  cellwise  velocity  components  are 
now  computed  using  equations  (4)  and  (5), 

Step  1,8  The  tentative  new  specific  internal  energy  is  computed 
using  equation  (6), 

Each  of  these  computational  steps  are  made  for  every  cell  in 
the  mesh  before  the  subsequent  step  is  made  for  any  one  cell. 

This  is  r  ssary  since  results  of  the  step  for  cell  A  are  usually 
required  t  compute  the  next  step  in  the  cells  neighboring  A, 

(b)  Phase  II  Calculations 

At  the  end  of  Phase  I  the  tentative  values,  o', ”v  and  I  are  known 
for  each  cell.  The  mass  particles  are  now  moved  and  this  effect 
is  taken  into  account  in  recomputing  these  qimntitics  in  the  steps 
depicted  in  Table  11  and  detailed  below. 
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Step  2, 1  This  is  a  preparatory  step  in  which  the  tentative 
energy  ana  momentum  are  totaled  over  all  the  particles  in  the  cell 
at  the  end  of  Phase  I: 


(17) 

il 

(18) 

11 

(19) 

II 

) 


Step  2,2  None  of  the  cellwise  quantities  are  changed  in  this 
computation;  the  mass  particles  are  moved.  If  r^,  are  the 
coordinates  of  a  mass  particle  at  the  end  of  the  n-th  time  cycle, 
its  new  coordinates  are  given  by 

r '  =  r  +  u  5t 
u  u  u 


(20) 


z  ' 
u 


Z  +  V  St  , 
I'  ^ 


where  and  are  computed  by  the  "area  averaging"  formulas  pre¬ 
sented  in  Appendix  B,  For  each  mass  particle  the  corresponding 

coordinates  r'.  z'  replace  (r  ,  z  )  in  the  machine  memory. 

1/  u  u  u 


step  2.3  The  total  mass  of  the  particles  in  each  cell  is  recom¬ 
puted  to  adjust  for  net  loss  or  gain  in  the  mass  present  in  the  cell 
after  the  motion  of  the  particles.  New  values  for  the  new  cellwise 
total  energy  and  new  values  for  the  components  of  the  total  cellwise 
momentum  are  also  computed. 


(21) 
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(22) 
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(24) 


zf  =  Z'  -  "i  S  %  +  E  f^N  E  N  ) 
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Here  "out",  means  that  the  summiition  is  taken  over  all  particles  leav¬ 
ing  cell  (i),  N  means  the  summation  is  taken  over  all  the  eight 
nearest  neighbors  of  cell^^^  and  the  "in"  means  the  mass  summation 
is  taken  over  all  particles 'entering  cell  from  a  given  neighboring 

cell  N,  Only  the  eight  nearest  neighbors  are  considered  since  St  is 
restricted  to  satisfy  condition  (B-1), 

Step  2.4  The  cellwise  density  and  the  velocity  components  at  the 


end  of 

the  (n  +  1)  -St 

(25) 

Mi' 

J  -  >  . 

p[ 

2  ff  Sr  Sz  ri 

1 

(26) 
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1  l/  1 
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Step  2,5  The  cellwise  specific  internal  energy  at  the  end  of  (n  +  1) 
-St  time  cycle  is  computed: 


(28) 


r-s-lto'.Or)'] 
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Step  2 . 6  The  pressure  at  the  end  of  the  (n  +  1)  -st  time  cycle 
is  computed: 

(29)  f  (p|  -  'I  ) 


At  the  end  of  Step  2.6  the  information  in  the  machine  memory 
is  in  such  a  form  that  the  computations  for  the  (n  +  2)  -nd  time 
cycle  may  begin  immediately.  However,  first  making  some  sub¬ 
sidiary  calculations  is  helpful. 

(c)  Phase  III  Calculations 

Before  proceeding  to  the  next  time  cycle,  computation  of  the 
total  energy  and  the  components  of  total  momentum  is  useful.  Thus, 
in  the  present  problem,  the  total  momentum  components  should  be 
rigorously  conserved,  but  the  total  energy  of  the  system  should  be 
a  monotone  decreasing  function  of  the  time.  Set 


R  =  S  hi 

(30)  System  ' 
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(32) 
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where  the  summation  is  extended  over  all  cells  containing  mass 
particles  disturbed  during  the  motion.  The  restrictions  on  the 
momentum  and  energy  then  become 


(33) 


R'  ^  R 
tot  ‘o‘- 
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(34) 


Z  =  Z 

tot  tot. 


Ii2 1 

—  H  I,.p  V 
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(35) 
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t.  =  F^tot.  =-".l!2l.p,v„2 


CONCLUSIONS 


In  applying  equations  (4),  (5)  and  (6)  in  Step  1.  1  the  values  of 
velocity  and  specific  internal  energy  must  be  known  in  neighboring 
cells  as  well  as  in  the  cell  in  question.  Specifically,  to  calculate 
'u!,  and  "v|  the  pressure  and  velocity  components  at  the  cell  centers 
and  botmdaries  depicted  in  Figure  3  mast  be  known.  The  usual 
computational  procedure  therefore  does  not  apply  near  a  free  sur¬ 
face  nor  in  the  vicinity  of  the  axis  of  symmetry.  Another  neces¬ 
sary  step  is  to  provide  for  such  contingencies  as  empty  cells  since 
the  usual  scheme  also  fails  for  that  situation. 

Current  effort  is  directed  toward  devising  means  for  handling 
these  special  situations.  It  is  too  early  to  report  on  these  tech¬ 
niques  in  detail;  it  should  be  possible  to  do  so  in  the  next  Quarterly 
Report,  Before  the  actual  programming  for  the  IBM  7090  is  com¬ 
menced  we  intend  to  visit  Dr,  Harlow  at  Los  Alamos  to  discuss  the 
scheme  in  detail. 
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APPENDIX  A 


Let  p,  p.ir  =  (u,  o,  v)  and  I  denote  the  density,  pressure,  velocity  and 
specific  internal  energy  respectively.  Then  the  axisymmetric  formu¬ 
lation  of  the  visco-plastic  equations  may  be  written  in  the  form 
(Ref.  1): 


(A-1) 


dp  dp  dp 

T+  “  T" 

dt  dr  dz 


+  p  div  u  =  0 


P 


+  u 


ilu  du 

dt  dz 


d  I  du 

-  P  +  2  p  (D) 


dr 


dr 


(A -2) 


+  2  P (D)  — 
dr 


dz 


P 


+ 


dv  dv  \ 

“  IT 


dv  ' 

P  +  2p(D)  — 
dz  , 


(A -3) 


/  de  de  de 

Pl-^  +  u— -+v  — 
\  dt  dr  dz 


J_  d 

r  dr 


u  P  +  V  Q  +  p  (D)— 
dr 


(A -4) 


dz 


+  — {vP+uQ  +  p  (D) 


d 

dz 


{A-5)  P  =  f(p,l) 
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where 


P  =  —  p - fi  (D)  div  u 

3 


0  =  p  (D) 


fz  (D)  =  ^  + 


/  du  <9v  \ 
y  dr  J 


V'd^ 


(A -6) 


The  components  of  the  stress  tensor  are  given  by 
r 

=  P  +  2  ,x  (D)  ^ 

^  r  =P  +  2p(D)-l 


(A -7) 


r  =  Q 

T’7  ^ 


and  the  flow  criterion  is  given  in  terms  of 


(A -8) 


=  T  ^+2  fi  T 
o  ”o  o 


Vi)^ 
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APPENDIX  B 


Consider  a  particle  in  cell(i).  During  the  (n  +  1)  -st  time 
cycle  it  may  remain  in  that  cell  or  it  may  move  into  one  of  the 
eight  nearest  neighbor  cells.  No  other  possibility  exists  since  we 
restrict  St  to  satisfy  the  inequalities 


(B-1)  St  u  <  Sr 

'  '  max 


St  V 


max 


<  5z  . 


Imagine  a  rectangle  of  cell  size  located  about  each  particle,  the 
particle  being  at  the  center  (see  sketch).  The  effective  velocity 
for  moving  the  particle  is  taken  as  the  weighted  average  of  the 
four  cellwise  velocities  of  the  cell  that  the  superposed  rectangle 
overlaps.  The  weighings  are  proportional  to  the  overlap  areas. 


^  r 


1 

1 

1 

1  » 

1  (r^,  z 

(15 . 

; 

-  -i 

The  cells  which  the  superposed  rectangle  overlaps  vary  with 
the  quadrant  of  cell  within  which  the  particle  lies.  The  pos¬ 
sibilities  and  corresponding  effective  velocity  components  of  the 
particle  are  as  follows: 
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First  Quadrant 


!(i  -  1/2)  5r  =  <  i  5r 

(j  -  1)  Sz  =  <  (j  _  1/2)  Sz 


{B-2) 


(B-3) 


Second  Quadrant 


(i  —  1)  5r  =  <  (i  —  1/2)  5r 

(j  -  1)  5z  =  <  (j  -  1/2)  Sz 


(B-4) 


U  = 


^  |[(i  -  1/2)  Sz  -  z J  [r^  -  (i  -  3/2)  Sr]  ui  -  '  +  [(i  -  1/2)  Sr  -  r J  uj  I 
[z^  -  (j  -  3/2)  Sz]  ^  [(i  -  1/2)  Sr  -  rj  u|  _  j  +  [r^  -  (i  -  3/2)  Sr]  u|  ^  | 
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t  "i  “ '  ^  ^  -  v]  "i  - 1) 

+  1^*1,  -  0  -  3/2)  j^(i  -  1/2)  ar  -  v|  _  J  +  -  (i  -  3/2)  arj  v|^  | 


Third  Quadrant 


(i  —  1)  ar  =  r^  <  (i  —  1/2)  ar 
(j  - 1/2)  az  =  <  j  az 


(B-6) 


"  ai^  I  ^  ~  ~  ~  [^' "  "’"*']  -  1  ^ 

+  -  (j  -  1/2)  azj  ^  j^(i  -  1/2)  ar  -  r  J  u|  ^  j  +  |^r^  -  (i  -  3/2)  arJ  uj  ■*■ 
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Fourth  Quadrant 


(i  —  1/2)  Sr=  <  i  St 
(j  —  1/2)  Sz  =  <  j  Sz 


(B-8) 


=  —^1  |^(j  +  1/2)  Sz  -  zj  ^  -  (i  -  1/2)  5rj  u|  ^  j  +  |^(i  +  1/2)  Sr  -  rj 

+  j^z^  -  (j  —  1/2)  Sz  ^  “  vj  “i  ^  +  j^V  ”  J  “i  +  1 


(B-9) 


[^(j  +  1/2)  Sz  -  zj  -  (i  -  1/2)  Sr]  v|  ^  j  +  (i  +  1/2)  Sr  -  r^ 

+  -  (j  -  1/2)  Sz]  ^  ^(i  +  1/2)  Sr  -  r^,]  v|  ■*■  *+  |^  r^^  -  (i  -  1/2)  Sr]  vl  J 


TABLE  I 
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|— H  — I 


Figure  1,  Illustration  of  projectile -tar get 
configuration  just  before  impact. 
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Figure  2.  Space  mesh  cells  of  PIC  method  superposed  over 
the  projectile-target  configuration  at  time  t  ■  o. 
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1 


Figure  3,  Layout  depicting  the  quantities  required  to  compute  tentative 

values  of  v*  and  7’. 

i  ’  i  i 
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